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Abstract. We present N-body simulations for generic non-Gaussian initial conditions with 
the aim of exploring and modelling the scale-dependent halo bias. This effect is evident 
O on very large scales requiring large simulation boxes. In addition, the previously avail- 

able prescription to implement generic non-Gaussian initial conditions has been improved 
Ci to keep under control higher-order terms which were spoiling the power spectrum on large 

scales. We pay particular attention to the differences between physical, inflation-motivated 
CN primordial bispectra and their factorizable templates, and to the operational definition of 

the non-Gaussian halo bias (which has both a scale-dependent and an approximately scale- 
rs) independent contributions). We find that analytic predictions for both the non-Gaussian halo 
mass function and halo bias work well once a fudge factor (which was introduced before but 
still lacks convincing physical explanation) is calibrated on simulations. The halo bias remains 
therefore an extremely promising tool to probe primordial non-Gaussianity and thus to give 
insights into the physical mechanism that generated the primordial perturbations. The sim- 
^-H ulation outputs and tables of the analytic predictions will be made publicly available via the 

J> non-Gaussian comparison project web site http://icc.ub.edu/~liciaverde/NGSCP.html. 
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1 Introduction 

Inflation, the standard theory for the origin of the primordial density fluctuations, predicts 
that the primordial density field can be described, to a good approximation, as a Gaussian 
random field. This nearly Gaussian nature of the primordial fluctuations is well established 
observationally (for current constraints see, e.g., [1]). However, small deviations from Gaus- 
sianity are predicted by inflation. The level and form of the non-Gaussianity depend on 
the inflationary model. The simplest inflationary model, the standard single-field slow-roll 
inflation, produces non-Gaussianities which are unmeasurably small [2, 3] (for a thorough 
review on inflationary non-Gaussianity see [4, 5]). 

Other inflationary models [6-12], which violate one or more assumptions of the standard 
single-field slow-roll inflation, may potentially cause much larger non-Gaussianities. Hence, 
measurements of deviations from Gaussianity offer a unique window to distinguish between 
different inflationary models or to rule out classes of models. 

Departures from Gaussianity are described at leading order by the bispectrum of the 
primordial fluctuations in the gravitational potential, B^{ki,k2,h^)- The bispectrum is the 
Fourier transform of the 3-point-correlation function and is identical to zero for a Gaussian 
field. The functional form of the bispectrum specifies the form of the non-Gaussianity, while 
its amplitude, parametrized by the parameter /nLj gives the level of non-Gaussianity. 

Since the Universe is statistically homogeneous and isotropic, the wave vectors ki , k2 ? 
and ks have to form a triangle. The shape of the triangle, for which the bispectrum peaks, 
is conventionally used to characterize the form of the non-Gaussianity. The local type of 
non-Gaussianity, ^{x) = ^^{x) + /nl [^^{x)"^ — (<I>"-^(x)^)) , where ^^{x) is gravitational 
potential in the Gaussian case, yields a bispectrum that peaks for squeezed triangles (e.g., 
ki <^k2 — k^). Potentially large non-Gaussianities of the local type are generally produced by 
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inflationary models with multiple fields, where non-Gaussianity is created by non-linearities 
which develop outside the horizon [8]. 

Other inflationary models which involve higher-derivative operators (see e.g. [13-15]) 
give rise to a bispectrum that has its maximum for equilateral triangles {ki = k2 = k^). 

Models in which the initial state of the inflaton is different from the Bunch-Davies 
vacuum [9, 11, 16, 17], can produce non-Gaussianity of the enfolded shape, i.e. the bispectrum 
is dominated by flattened or enfolded triangle configurations {ki ~ A;2 + k^). 

Finally, there is the so-called orthogonal shape of non-Gaussianity, which cannot be 
easily described in terms of a triangle configuration [18]. This shape is — with respect to a 
specific scalar product (see [8] for an explicit definition) — orthogonal to the equilateral and 
local shape. 

The cosmic microwave background (CMB) offers the cleanest probe of non-Gaussianity, 
since the temperature perturbations are still very small and well described by linear theory. If 
a shape of the primordial bispectrum is assumed^ , constraints on the corresponding /nl value 
can be derived from the observed fluctuations in the CMB. The current constraints for the 
most commonly used shapes are -10 < < 74, -214 < f^f'^ < 266, -410 < < 6 

(95% CL) [1], and f^""^ = 114 ± 72 (68% CL) [20] . The Planck mission [21] is expected to 
measure /nl with an accuracy of A/nl — 5. 

Primordial non-Gaussianity leaves also observable signatures in the large-scale structure 
of the Universe (for recent reviews, see [22, 23]). The bispectrum of the density field measured 
through observable tracers like galaxies, however, is dominated by non-linear gravitational 
evolution and biasing, this makes the measurement of the primordial contribution difficult 
[24-26] (see, however, [27]). 

Another probe of non-Gaussianity from large-scale structure observations, is the abun- 
dance of galaxy clusters or large voids in the Universe. Both of these objects originate from 
the tails of the nearly Gaussian initial density distribution and are therefore sensitive to pri- 
mordial non-Gaussainity [28, 29]. The signal depends on the skewness of the initial density 
distribution, which is given by an integral over the bispectrum, and is therefore less sensitive 
to the shape of the non-Gaussianity. This probe measures the primordial non-Gaussianity on 
scales of the order of lOMpc, which is smaller than the scales accessible by CMB measure- 
ments. Intriguingly, several groups [30-32] derived estimates of /nl from the observations of 
very massive galaxy clusters, which are as large as /nl ~ 400 (see, however, [33]). 

The clustering of halos on large scales offers another possibility to probe non-Gaussianity 
[34, 35]. Non-Gaussianity couples small and large scales. The density peaks on small scales, 
which are the seeds for halo formation, are modulated by the large scale modes of the grav- 
itational potential. This induces a scale-dependent halo bias on large scales, which can be 
very different from the scale-independent halo bias in the Gaussian case. The amplitude 
of this effect scales approximately linearly in /nL; while the scale dependence is set by the 
shape of non-Gaussianity. For example, the local type non-Gaussianity generates a halo bias 
which scales as k~'^, whereas the halo bias of equilateral types of non-Gaussianity is expected 
to have only a very mild scale dependence [36]. This technique was already used on ob- 
servational data of galaxy surveys and quasar catalogues [37-39]. The derived constraints, 
-27 < /^£^^ < 70 (95% CL) [37], are competitive with the CMB measurements. Fisher 
matrix forecasts for future large-scale structure surveys predict that measurements of the 
non-Gaussian bias will constrain /^£^' with an uncertainty of A/nl ~ 1 [40-43]. 

^For an approach without assuming a shape for the bispectrum, see [19]. 
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All of the three aforementioned large-scale structure probes of non-Gaussianity are 
highly affected by non-linear gravitational evolution. The standard tool for taking non- 
linearities into account are N-body simulations. However, until recently, only the local type 
of non-Gaussianity was simulated with N-body simulations and reasonable agreement with 
analytic predictions was found [34, 44-47]. 

In [48], we presented a prescription for setting up initial conditions for N-body simu- 
lations with non-local non-Gaussianities. Using this method we were able to simulate the 
non-linear power spectrum and the cluster mass function for non-local shapes, and again 
found consistency with theoretical expectations. The technique of [48] was, however, not 
suitable to study the effect of non-Gaussianity on the halo bias, since — depending on the 
shape of the bispectrum — it introduced artificial power on large scales, which interferes with 
the expected signal of the non-Gaussian halo bias. 

In this paper, we modify our ansatz used for the generation of non-Gaussian initial 
conditions in a way that these spurious contributions to the initial power spectrum are always 
kept under control. This enables us to set up and run non-Gaussian simulations with large 
volumes, which we then use to study — for the first time by means of N-body simulations — 
the non-Gaussian halo bias for different types of non-Gaussianity. As the simulated volume 
is larger than the one we were able to simulate in [48] and hence the statistical error on the 
number of massive halos is smaller, we use these simulations also to revisit the non-Gaussian 
cluster mass function. 

The outline of the paper is the following. First, we review the non-Gaussian halo bias 
and its theoretical modelling in Sec. 2. In Sec. 3, we compare, in the context of the non- 
Gaussian halo bias, the bispectra of inflationary models with the commonly used separable 
templates which approximate these bispectra. In Sec. 4, we review the method of setting up 
non-Gaussian initial conditions for a given bispectrum and then present the modification of 
our ansatz. The practical implementation and numerical settings of our simulation suite are 
given in Sec. 5. We present and discuss our results on the non-Gaussian mass function and 
halo bias in Sec. 6. Finally, we conclude in Sec. 7. 

2 Non-Gaussian halo bias 

The halo bias describes the relation between the halo density contrast, Si^, and the un- 
derlying dark matter density contrast, (5ni- In Fourier space we can write the relation as 
^hik) = b{k)5ra{k) + n{k), where the random variable n{k) models the stochasticity coming 
from shot noise and possibly from physical halo formation processes. For Gaussian initial 
conditions the bias is scale-independent on large scales, bM,z{k) = bf-M z ^ ^ ^-^ [49-51]. 
Analytic predictions and results of N-body simulations show that the linear^ Gaussian bias, 
^i^Mz' increases with redshift z and halo mass M. In addition the halo bias may depend on 
environment and mass accretion history [52, 53]. 

Primordial non-Gaussianity of the local type induces a scale dependence of the halo 
bias on large scales, which scales as ~ This particular scaling was first found in N-body 
simulations by [34], who also gave a simple analytic understanding how this effect arises. 
Later [35, 37] rederived the analytic prediction in a more rigorous way using two different 
approaches, the high peak limit and the peak-background split, respectively. In Ref. [54], the 

■^The linear bias, bi, is called linear as it corresponds to the linear coefficient in the Taylor expansion of 
the halo density field: Sh{x) — bo + bi5m{x) + fe25m(a^)/2 + . . . , where both density fields are smoothed on a 
sufficiently large scale. 
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same scale dependence of the non-Gaussian halo bias was found by computing the ellipsoidal 
collapse threshold in the presence of small non-Gaussianities. The analytic predictions were 
tested with further and larger N-body simulations [44-47]. Reasonable agreement between 
theory and simulations was found if small corrections, which we discuss in detail below, 
are taken into account. In a number of subsequent papers further theoretical modelling 
was developed by applying (univariate and multivariate) local biasing in the framework of 
perturbation theory [27, 55-59]. In summary the non-Gaussian bias can be modelled to 
leading order as 

''(2.1) 

where A6i;M,2(/nl) is the difference in the linear bias, which is non-zero, since a non- vanishing 
/nl changes the halo number density and hence the linear halo bias for a fixed halo mass M. 
Mainly, this effect leads to a scale-independent shift relative to the Gaussian (/nl = 0) linear 
bias [44, 59, 60]. Additionally, it makes the scale-dependent second term in above equation 
slightly non-linear in /nl [59, 61]. In practice, when fitting observational data, the linear 
bias is estimated from the data itself. 

The amplitude of the non-Gaussian contribution is given by /nLj the linear Lagrangian 
bias [&1;A//,z(/nl) — 1], and the redshift-dependent collapse threshold q6c/D{z), where 5z = 
1.686 is the spherical collapse threshold in an Einstein-de Sitter universe, D{z) is the growth 
function normalized to (1 -|- z)~^ at high redshift, and q can either be regarded as a fudge 
factor introduced to improve the fitting to N-body simulations [46] or as a correction of 
the spherical collapse threshold [54, 62]. The scale dependence of the non-Gaussian bias is 
hidden in A^Af(^) and FM{k)- M.M{k) connects the density with the primordial potential. 
The linear matter density contrast smoothed on a mass scale M is related to the primordial 
potential by the Poisson equation, 

2T{k)k'^ 

Here, and Hq are the present day matter fraction and the Hubble constant, respectively. 
The transfer function T(k) models the physics until the end of recombination and is nor- 
malised to unity on large scales. The function W^m(^) is the Fourier transform of the top-hat 
filter with a smoothing scale given by the mass M. As W^f{k — )• 0) = 1, we see that AiM^k) 
scales as k'^ on large scales. 

The "form factor", Tjviik)-, depends on the shape of the non-Gaussianity through the 
bispectrum of the primordial potential B^{ki, ^02,^3), [35]: 

J'Mik) = ^ ,1 ^ I dk'k'^Muik') C d^MM{\k + \^\)^^^^^4^%^^ ' (2-3) 

where ^ is defined by k • k' = kk' ^i, P^{k) is the primordial power spectrum, and a^.^ 
denotes the variance of linear density fluctuations on the mass scale M at redshiht z = 0. 
Considering the limit of /c — ?■ 0, we see that on large scales the form factor depends primarily 
on the squeezed limit of the bispectrum. For the local type of non-Gaussianity, the "form 
factor" becomes scale and mass independent on large scales, FM{k ^ 1) = 2. Hence, we 
obtain the aforementioned k~'^ scale dependence of the halo bias in the case of local non- 
Gaussianity. As we discuss below in Sec. 3, other shapes of non-Gaussianity lead to different 
and in general weaker scale dependences. 



5m;M,.(k) = --^-^WM{k)D{z)<^{\L) = MM{k)D{z)^{k) . (2.2) 
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One higher-order correction to A6, which was introduced in [44], is a relatively small 
contribution coming from the fractional difference in the Gaussian and non-Gaussian non- 
linear matter power spectrum. This effect becomes larger at smaller scales but remains below 
a few percent on all scales for observationally allowed values of /nl- 

Some additional effects are also missed in Eq. (2.1): The impact of the mass accretion 
history of the halos — in particular recent major mergers — on the non-Gaussian halo bias 
was pointed out in [37], and further investigated and tested with N-body simulations in [63]. 
The formation redshift of the halos, zj, changes the amplitude of the non-Gaussian correction 

by bf.M -^ifNL) - 1 — > bf.M^zifNL) - K^f) , where iJ,{zf) varies roughly between -2 and 
+2 from very young to very old halos [63]. For recent major mergers, /i takes the value 
-1 - l/6c ~ -1.6 [37, 63]. 

In addition, higher-order effects in the biasing description (e.g., terms of the order /^^ 
[55, 59, 64]) and in the primordial non-Gaussianity (e.g., non-Gaussianity at the cubic order 
quantified by ^nl [59, 64]) are neglected in Eq. (2.1). 

Finally, general relativistic effects are not taken into account, which might become 
important on very large scales [36, 65-69]. 

So far we discussed the scale-dependent halo bias caused by the local type of non- 
Gaussianity with a scale- independent /nl- The effect of a scale-dependent /nl were inves- 
tigated analytically and by means of N-body simulations in [70, 71]. In this paper we are 
mainly interested in the scale-dependent halo bias of non-local types of non-Gaussianity. The 
prescription presented above is sufficiently general to model the effects of other shapes of non- 
Gaussianity. As discussed above, the shape dependence is hidden in the form factor FM{k)- 
The expression for the form factor, see Eq. (2.3), was derived by [35] in the framework of 
local biasing of density peaks and was evaluated for different types on non-Gaussianity in 
[36]. Recently, [72] generalized the peak-background split approach to non-local types of 
non-Gaussianities and derived a very similar but not identical expression for the form factor. 
In the next section we discuss the analytic predictions of the non-Gaussian halo bias for 
non-local shapes using the high peak formulation of [35] . 



3 Physical shapes vs. templates 

The functional form of the primordial bispectrum, B^{ki,k2-,k^)-, predicted by inflationary 
models can be quite complex. In particular, the bispectrum is usually not separable, i.e. it 
cannot be written as a finite sum of terms which are factorizable as a product of functions of 
ki, k2 and k^. Separability is important for computationally efficient analysis and simulation 
of both CMB maps [13, 19] and large-scale structure [48, 73]. Hence, most physical shapes 
have been approximated by so-called templates, which are separable and constructed to 
effectively maximize the correlation between the physical shape and the template across all 
triangle configurations. In practice a so-called "cosine" shape correlator is used [8], which 
absolute value is always < 1 with 1 indicating perfect correlation. This measure of similarity 
is useful for constraints on the bispectrum coming from CMB or large-scale structure analysis, 
as it probes the entire range of possible triangle configurations. However, these templates 
can be misleading for predictions of the non-Gaussian halo bias, as in this case the shape 
is probed mostly in the squeezed limit. In this limit the templates do not necessarily need 
to be good approximations of the physical shapes to yield a cosine close to unity [20]. For 
the non-Gaussian halo bias, however, the correct scaling of the templates in this regime is 
crucial for sensible predictions. Let us consider several templates in this regard. We start 
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Figure 1. Scale-dependent part of the non-Gaussian lialo bias, TM{k)/MM{k), for the equilateral 
types of non-Gaussianity. Solid and dashed lines correspond to a spectral index of Us — 0.95 and 
rig = 1, respectively. In both cases, the amplitude of the power spectrum is set by keeping the mass 
variance fixed, as — 0.7913. 

with the equilateral'^ template, which was introduced in [74] as a separable approximation 
to the shape functions of DBI inflation [14], inflation with higher-derivative terms [13], and 
ghost inflation [15]: 

Bfih, k2, h) = 6f^i - 2F^ + F^) , (3.1) 

with^ 

f1°^ = P^{ki)P^{k2) + P^{k2)P^{k3) + P^{ki)P^{k3) (3.2) 
= [P^{ki)P^{k2)P^{k3)f^ (3.3) 
= Pl^\ki)Pl^\k2)P^{k3) + 5perm. (3.4) 

and the dependence of F on the three wavenumbers k is understood. The primordial power 
spectrum is specified by the amplitude at /cq and the spectral index n^: 

P^{k) = Ao {k/koT'-^k-^ . (3.5) 

The equilateral template was also adopted in other models. More generally, the form 
of the bispectrum is interesting because of its direct connection to the (self) interactions 

^This template peaks when the wavenumbers k form a equilateral triangle. 

^The original expression of the equilateral template was proposed for a scale-invariant power spectrum 
and given in terms of the wavenumbers k. Here we generalize it to non scale-invariant power spectra by 
substituting k with P~^^^{k) with proper normalization. We use the same substitution ansatz in the case of 
other bispectra. 
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of the inflaton. In the effective field theory of inflation [75], there is a direct connection 
between the coefficients of operators of the Lagrangian for the inflationary fluctuations and 
the shapes of primordial non-Gaussianity. In [18], it was shown that the bispectrum of single- 
field inflation with an approximate shift symmetry protecting the Goldstone boson is a linear 

combination of two different shapes. The two shapes, B'^^^'''^^ and -BJ^, which correspond to 
the interaction operators 7r((?j7r)2 and vr^ (see [18] for details), are both close to the equilateral 

template. Especially, the cosine of * with the equilateral template is very close to one. 
For this reason, we denote this shape also by "SSZ equilateral" 

We now turn to the question how well the equilateral template captures the behaviour 
of these physical shapes in the squeezed limit. Choosing fci — )• and using the following 
notation, k = ki, k' = k.2, and = |k + k'| = \/k'^ + fc'^ + 2k'kn, we obtain for the 
equilateral template 



^eql _ 12 f'^'^^ 

^ squeezed •' NL 



pi^\k)Pi/\k') - (^"^y t,'p^{k)k'p^{k')k'-' 



(3.6) 



NL 



' kh'l V' / n,-4 y ^fkt 



k-^k'-^ (3.7) 



= 12Alf^t{l-^^')k-'k'-'^k-\ (3.8) 

where in the last line a scale-invariant power spectrum was assumed, i.e. Ug = 1. We 
find a similar scaling of the bispectrum in the squeezed limit also for the physical shapes 
generated by the aforementioned inflationary models. However, the exact scaling and its 
amplitude differ slightly from case to case. In the case of a scale-invariant power spectrum, 
the scaling reduces for all models to k^^ with a model-dependent amplitude. For example 
for the bispectrum of DBI inflation (see [8] for an explicit expression) we get 

nDBI _ 99 /l2fcql f , 5 o] j,-lj,'~5 (o q\ 

-D* squeezed " y-^O/NL j K. ^ , ['i-y) 

while the scaling of the bispectrum corresponding to n^diTr)'^ [18] is 

,SSZ eql. _ 216 2 .eql 5 o ^ , -1, /-5 



Br ''^'•squeezed = -^^tf^l " ^^^J ^'^^'"^ • (3-10) 

Here, we used the equilateral normalization convention for all shapes, i.e. B^{ko, ko, ko) = 

K 7121,-6 ^eql 

In Fig. 1, the scale-dependent part of the non-Gaussian halo bias, FM{k) / Ai^ik), is 
shown for different bispectra of the equilateral type. The dashed lines depict J^Mik)/M.M{k) 
for a scale-invariant power spectrum. Note that on very large scales the amplitude of the 
effect for the equilateral template and the physical models, namely DBI and SSZ eql., is 
different by a factor of 1.5 and 1.256, respectively, in agreement with the above expressions 
for the squeezed limit. The solid lines show J^a/ (/c) /M M{k) computed with a power spectrum 
with Ug = 0.95. The upturn and the wiggles on smaller scales stem from the transfer function 
included in M.M{k). Note, in contrast to the local type of non-Gaussianity, TM{k) / M.M{k) 
for the equilateral shapes is mass-dependent even on large scale. 



^Previously, in the framework of general single field infiation, [9] obtained two independent shapes which 
span the same space as the two shapes of [18]. 



-7- 



In summary, predicting the non-Gaussian halo bias with the equilateral template instead 
of the physical shapes results in approximately the correct scaling of the effect but a slightly 
smaller amplitude. The amplitude derived from the physical shapes of DBI and SSZ eql. is 
larger by 50% and 26%, respectively. 

For other templates, as we discuss now, the situation becomes more problematic. The 
orthogonal^ template [18] 

Bl''\ki, k2, h) = df^f" {-3F'°^ - 8F^ + 3F«) , (3.11) 

and the enfolded^ template [16] 

scale both as ~ k~'^ in the squeezed limit, while the corresponding physical shapes scale 
as ~ k~^ and ~ k~'^, respectively. In Fig. 2 we show again the scale-dependent part of 
the non-Gaussian bias, FM{k) / M.M{k), computed from these commonly used templates and 
the corresponding physical shapes. For comparison, the scale dependence obtained from the 
local type of non-Gaussianity, B^^^ = 2f^^F^°^, and from the equilateral template are also 
shown. Here, we used in all cases a scale-invariant power spectrum. The solid and dashed 
lines correspond to halos of mass 3 x IO^^'Mq/Zi and 3 x IO^^Mq/Zi, respectively. The scale 
dependence of the halo bias caused by modified-initial-state type of non-Gaussianity [16] 
is almost degenerate with the one caused by the local type of non-Gaussianity, i.e. both 
predictions are very similar when is rescaled by 1/8. However, on scales of ~ 0.01/i/Mpc 
the different halo mass dependence lifts this degeneracy to some extent. The orthogonal shape 
of [18] (SSZ orth.) — having the same scaling in the squeezed limit as the equilateral shapes — 
does not lead to a scale-dependent halo bias on large scales and is almost indistinguishable 
from the equilateral type. Only on smaller scales, the "orthogonal" halo bias has a slightly 
different scale dependence than the equilateral one. 

In conclusion, in respect of the halo bias, the commonly used enfolded and orthogonal 
templates are not adequate to capture the behaviour of the underlying physical models. The 
physical shapes scale as ~ k~^ or as ~ k~^ in the squeezed limit, while the scaling of the 
enfolded and orthogonal templates is ~ /c~^. 

In [19] a mode expansion of the bispectrum in separable functions was introduced. This 
expansion can be applied to arbitrary non-separable bispectra. In many cases already a 
small number of modes is sufficient to yield a large cosine with the original bispectrum. 
While this method of constructing separable templates for arbitrary bispectra is useful for 
efficient measurements of the bispectrum from the CMB or the large-scale structure, it fails 
to reproduce the correct scaling of the bispectrum in the squeezed limit. By construction of 
the mode expansion, all templates derived in this way have a scaling of ~ k~'^ in the squeezed 
limit, irrespectively of the scaling of the original bispectrum. 

New separable templates which are not only overall good approximations to the original 
bispectrum but also faithfully capture the scaling in the squeezed limit can most probably 
be constructed (e.g [20]). However, we will see in the next section that separability of the 

®The orthogonal template was introduced as a separable approximation to the linear combination of the 
shapes B^'^''^^ and that is orthogonal (i.e. has a vanishing cosine) to B^^^^^^ 

^This template approximates the shape function produced by inflationary models with a modifications of 
the initial-state of the inflaton field and has its maximum for flattened triangles, i.e. fci = k2 + kg. 
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Figure 2. Scale-dependent part of the non-Gaussian halo bias, J-'M{k)/MM{k), for different infla- 
tionary models and their corresponding templates. The solid and dashed lines show the results for 
halos of mass 3 x 10^'* Mq/H and 3 x 10^^ Mq/H, respectively. 



bispectrum does not speed up our current implementation of generating the initial conditions 
for N-body simulation. 

As a concluding remark, we note that similar considerations apply for the skewness 
specified by the bispectrum. The skewness is the relevant parameter used for the predictions 
of the non-Gaussian halo mass function. Although having a large mean correlation over all 
triangle configurations, templates and corresponding physical shapes may lead to different 
values of the skewness. However, the differences in the skewness will be, in general, much 
smaller than the aforementioned differences in the scaling behaviour of the squeezed limit. 

4 Initial conditions 

In [48] we introduced a prescription how to generate non-Gaussian initial conditions for a 
given bispectrum, i.e. how we find a realization of the potential field which has the desired 
bispectrum 

(^ki^k^^-ka) = (27r)35D(ki + k2 + k3)S$(ki,k2,k3). (4.1) 

We shortly review our method. First, we decompose the potential field in Fourier space, ^>k, 
in a Gaussian and a non-Gaussian part 

$k = ^>k + ^'k^ • (4.2) 

The Gaussian field is statistically completely described by the power spectrum 

i^l^l) = (27r)35^(ki + k2)P^{h) . (4.3) 
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As the non-Gaussian contribution is small, the leading order of the bispectrum is given by 



(^ki^k.^-ka) = (^^1^,^^) + i'^K'^'^Z) + ■ (4-4) 

Using the following ansatz for <I>^'^ 

K2 K3 



= 6ihf I I' + "^i' ■ <^-^' 



^P<i.(A:2)P$(A;3) 



with <I>^ the complex conjugate of $k, one can easily show that 

^3«) = ^-{2T,fB^{ki,k2, A;3)5^(ki + k2 + ka) . (4.6) 

Finally, we recover at leading order Eq. (4.1) by virtue of the symmetry property of the 
bispectrum. 

One problem of the above ansatz is that the non-Gaussian contributions to the power 
spectrum can become very large on large scales. ^ Naively, one expects {(^^^(^^^) to 
be much smaller than the Gaussian contribution ($'^$*^). However in the case of certain 
bispectra, diverges more strongly for small wavenumbers than the Gaussian power 

spectrum reads 



spectrum. With the above ansatz for ^jT^ the non-Gaussian contribution to the power 



(^k )-^^^ ^^ + "1^ P<i.(fcO^<i.(|k + k'|)- 

In practice, there is a minimum wavenumber /cmin; the largest mode which fits in the simula- 
tion box, and a maximum wavenumber fcmax) which corresponds to the smallest wavelength 
still resolved by the used grid size. These wavenumbers provide large and small-scale cutoffs 
for the above integral. Nevertheless, this spurious contribution to the power spectrum can be 
large enough to affect substantially the results of the N-body simulations. In particular, sim- 
ulations with large box sizes, which are needed to predict the halo bias on large scales, suffer 
from these divergences. On large scales. A; ^ 1, the integral is dominated by the bispectrum 
in the squeezed limit. Using the scalings presented in Sec. 3, we see that the non-Gaussian 
contribution scales as f^^kP^{k), f^^k~^ P^{k), and f^i^k~^P^{k) for the equilateral, or- 
thogonal/enfolded templates, and the local type of non-Gaussianity, respectively. Hence, in 
the case of the equilateral type the contribution is small, whereas in the other cases it starts 
to dominate the power spectrum on large scales.^ 

For the local type of non-Gaussianity, we found in [48] that our ansatz becomes identical 
to the real space formulation <I)(x) = (^"^(x) + /nl [(<^*^(x))^ - ((<I>'^(x))^)] , which does not 
lead to the aforementioned divergences, when we use B^{ki, k2, k^) = 6/nl-P^(^2)-P$(^3) in 



On large scales, we require the non-Gaussian contributions to the power spectrum to be negligible (or 
renormalizable) , in order to be in agreement with CMB measurements. 

^In [73] the mode expansion technique of [19] was used to develop a computationally efficient method to 
generate initial conditions for N-body simulations for a given bispectrum. However, due to the k~'^ scaling 
in the squeezed limit of the expanded bispectrum, also in this case, non-Gaussian contributions to the power 
spectrum of the form f§ji^k~^ P^{k) arise on large scales. 
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Eq. (4.5) instead of B^^'^. Here we generalize this modification by writing it in the fohowing 
way 

Inserting this expression in Eq. (4.5) yields the modified ansatz for ^j^*^ (see also [72]) 
^NO_^f^3j^, B^ik,k',\l. + l.'\)<^t9<^o^^, 



2(27r)3 J P^{k)P^{k') + P^{k')P^{\k + k'l) + P^{k)P^{\k + k'| 

which by using Eq. (4.4) reproduces Eq. (4.1). 

Note that the integrand in the modified ansatz is in general not separable. This means 
that the integration cannot be written in a simple way as a sum of convolutions, which 
can be efficiently computed by Fast Fourier Transformations (FFTs). However, a possible 
avenue to still factorize our modified ansatz is the use of Schwinger parameters (see [76] 
for details). This is in contrast to our original ansatz, which is separable as long as the 
bispectrum is separable. However, the modified ansatz does not lead to contributions to the 
power spectrum, which diverge for A; — )• 0, since they are suppressed by P^ {k)-}^ 



«^0 = ^^^(k + q)/rf^fc' 



l.Dn, , „^ f Bl{k, k', |k + k'\)P^{k')P^i\k + k^l) 

[P$(A;)P$(fc') + P^{k')P^{\k + k'l) + P^{k)P^{\k + k'\)f ■ 

(4.10) 

Note that neither Eq. (4.5) nor Eq. (4.9) specify the trispectrum of the non-Gaussian 
field at the leading order. A possible inclusion of the trispectrum in the generation of the 
initial conditions was proposed in [73]. In this paper, however, we neglect contributions from 
higher-order correlations. 

Next, we can convert the potential into the linear density field with the help of 
the transfer function, which we compute with CAMB [77], and the Poisson equation: 

s^^ifmrn, ,4.11) 

3 ilmHo 

Finally, the density is sampled with particles. The positions and velocities of the particles are 
derived from the displacement field at the initial redshift using the Zel'dovich Approximation 
or second-order Lagrangian perturbation theory [78-80] . 



5 Simulations 



In this paper we are mainly interested in the scale-dependent halo bias of different types of 
non-Gaussianity. As this effect, in general, is larger on large scales, we need to simulate large 
volumes but still have sufficient mass and force resolution to resolve halos. Since we are not 
interested in halos of a particular mass, a reasonable trade-off between size and resolution is 
a box size of about 2Gpc/h and a particle load of one billion particles. This enables us to 
probe large scales {k ~ 0.003 /i/Mpc) and resolve halos with masses above lO^^Mo/Zi (~ 20 
particles). 

^"Nevertheless, without cutoffs the integral can still be divergent, e.g. in the local case this would lead to 
a renormalization of the amplitude. In practice, the cutoffs given by the box size and grid size render this 
contribution negligible. 
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The downside of our modified ansatz for the initial non-Gaussian part of the potential, 
given in Eq. (4.9), is the non-separability of the integrand. This precludes us from using 
FFTs to compute the integral. Instead, we perform the integration in Fourier space by direct 
summation, i.e. for each k mode the following sum has to be computed 



^NG 



k+k' 



2 ^P^{k)P^ik') + P^{k')P^{\k + k'l) + P^{k)P^{\k + k'l 



(5.1) 



where k' runs over all grid points in Fourier space. If Ng denotes the number of grid points 
in one dimension, then the computational cost of computing for all k modes scales 

as Ng. This scaling renders the usage of large grid sizes infeasible. Even for a moderately 
small grid size of 400, the generation of the non-Gaussian initial conditions takes two days 
on 256 cores of present-day CPUs. In order to still use a large particle load, we use two 
different grid sizes: One grid being as large as the particle grid, Np, which samples only the 
Gaussian contribution to the potential, and one smaller grid for the non-Gaussian part.^^ 
By doing this, we neglect non-Gaussianities on scales smaller than the grid spacing of the 
non-Gaussian grid. This not only changes the clustering on these scales, but may also affect 
the halo clustering on large scales, as the formation of the halos occurs on small scales. This 
effect will depend on the halo mass; halos which form on scales larger than the resolution of 
the non-Gaussian grid will, presumably, hardly be affected. The choice of the non-Gaussian 
grid size is therefore a trade-off between computational resources and the lowest halo mass, 
which is still unaffected by the spatially unresolved non-Gaussianity. 





k [h/Mpc] 



k [h/Mpc] 



Figure 3. Halo bias for different halo masses computed by the ratio of the halo-matter cross 
power spectrum to the matter power spectrum obtained from N-body simulations with local non- 
Gaussianities with /nl = 250. The solid and dashed lines correspond to local- 1024- 1024 and 
local-1024-400, respectively. The masses are given in units of Mq/H. 

To quantify this effect in more detail, we run two simulations of the local type of non- 
Gaussianity with /nl = 250. The settings of the simulations are identical except for the 
way we set up the initial conditions. In one case, we use a Gaussian grid of size 1024 and 
a non-Gaussian grid of size 400 (local-1024-400), in the other, both grids are of size 1024 
(local- 1024- 1024). This large non-Gaussian grid size of 1024 is possible because in the 

^^To compute the non-Gaussian part, we use Eq. (5.1) where we only sum over modes which he within the 
smaUer grid. 
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case of local non-Gaussianity the computation can be performed efficiently in real space, 




Figure 4. The ratio of the number density of halos found in the simulation local- 1024-400 and in 
the simulation local-1024-1024 as a function of halo mass. The error bars show the Poisson error. 
For clarity, the green circles are shifted slightly along the x-axis. 

We adopt a flat ACDM cosmology with cosmological parameters consistent with current 
observational constraints [1], namely ilm = 0-27 (the present-day matter fraction), ilf, = 0.047 
(present-day baryon fraction), Hq = 70km/sMpc~^ (the Hubble constant), Ug = 0.95 (the 
spectral index of the primordial power spectrum), and as = 0.7913 (the rms of the linear 
density fluctuations at z = in spheres of 8Mpc//i). 

The simulations are carried out with Gadget-2 [81] taking only the gravitational inter- 
action into account. The box size is 1875Mpc//i, and 1024'^ particles with a softening length 
of 40kpc//i are used to sample the density field. The simulations are started at an initial 
redshift of ^^mitiai = 49 using second-order Lagrangian perturbation theory to displace the 
particles from a regular grid. 

Halos are found in the particle outputs of the N-body simulations by the publicly avail- 
able halo finder AHF [82]. This halo finder defines halos as gravitationally bound objects 
with a spherical overdensity equal to the redshift-dependent virial overdensity. With the halo 
catalogues and particle snapshots at hand, we can compute the Fourier modes and 5k of 
the halo and matter density field, respectively. This is achieved by first assigning the halos or 
particles to a regular grid (512^) using the cloud-in-cell (CIC) scheme and then performing 
an FFT. 

In Fig. 3, we show the halo bias obtained from the two simulations local-1024-1024 
(solid lines) and local- 1024-400 (dashed lines) for different halo masses at redshift z = 
and z = 1. We compute the halo bias by the ratio of the halo-matter cross power spectrum 
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Table 1. N-body simulations. Ng denotes the size of the grid used for 
the non-Gaussian part of the potential. The size of the particle grid 
is identical to the size of the grid used for the Gaussian part of the 
potential and is given by Np. 



type oi non-Gaussianity 


JNL 


A 7 

^9 


Np 


7f realizations 


local 


250 


400 


1024 


1 


local 


250 


1024 


1024 


2 


local 


60 


1024 


1024 


3 


equilateral template 


1000 


400 


1024 


2 


equilateral template'^ 


1000 


1024 


1024 


1 


orthogonal template 


-1000 


400 


1024 


2 


orthogonal template 


-250 


400 


1024 


3 


Gaussian 






1024 


3 



"In this case, we use ansatz Eq. (4.5) for the non-Gaussian . The 



convolution is computed with an FFT in real space. 



to the matter power spectrum 

with P^{k) = (|JkP) and Phm(fc) = (Re(5jj(5^)), where the average is taken over all modes 
which lie in a shell of radius k and thickness AA; = 0.003 /i/Mpc. As expected, the clustering 
of halos with a mass below ~ 3 x IO^^Mq/Zi is affected by the smaller grid size used for 
the non-Gaussian <I>^^. The grid spacing of the non-Gaussian field is 4.7Mpc//i. This 
corresponds approximately to the Lagrangian radius of halos of mass 3 x lO^^Mo/Zi. The 
small difference in the halo clustering for more massive halos is consistent with shot noise. 
At higher redshift, the noise increases due the smaller number of halos per mass bin. 

The other quantity of interest, which is also affected by the smaller non-Gaussian grid, 
is the halo number density. In Fig. 4, we show the ratio of the mass functions derived from 
the two simulations local-1024-400 and local-1024-1024. We find that the smaller non- 
Gaussian grid size decreases the halo number density at the low-mass end by a few percent. 
For halos more massive than ~ 4 x 10^^ Mq/Zi, the halo number densities of both simulations 
agree within the errors. 

In App. A, we compare the results of two N-body simulations of the equilateral type, 
for which we used the two different formulations for the non-Gaussian part of the potential, 
Eq. (4.5) and Eq. (4.9). This comparison gives us another estimate of the minimum halo 
mass, ~ 5 X 10^^ Mq//i, above which the results are not affected by the low spatial resolution 
of the non-Gaussianity. 

We conclude that, for halos more massive than 5x10^3 Mq/H, even the simulations which 
are set up with a small non-Gaussian grid {Ng = 400) model correctly the characteristics 
of non-Gaussianity, in particular the halo bias and the halo number density. As a rule of 
thumb, the non-Gaussian grid spacing should be smaller than the Lagrangian radius of the 
lowest-mass halo of interest. 

For the main study of the non-Gaussian halo bias, we perform a suite of simulations 
of different shapes of non-Gaussianity and different /nl values, which are summarized in 



- 14 - 



Tab. 1. The cosmology and simulation settings, like box size, particle load, etc., are the same 
as given above. 

Although the orthogonal template is not a good approximation of the physical shape 
(SSZ orth.) in the context of halo bias (see Sec. 3) and although our ansatz is applicable 
to non-separable bispectra, we still use this template for some of our simulations. This has 
the following reasons. First, the corresponding physical shape (SSZ orth.) causes a halo bias 
which is very similar to the equilateral type and, hence, hard to measure. Second, the scale 
dependence of the non-Gaussian halo bias of the orthogonal template lies in between the 
local and the equilateral one, which makes it an interesting toy model to test the analytic 
predictions of the non-Gaussian halo bias against N-body simulations. Lastly, other groups 
(e.g. [83]) are running N-body simulations of this type, too, which makes a future comparison 
easier. 

The chosen values for /nl are partly larger than allowed by current constraints coming 
from CMB data [1]. Here, however, we want primarily to compare theoretical predictions 
with the results of N-body simulations. To do this accurately, we need high signal-to- noise 
and therefore a relatively large value of /nl- This is in particular the case for the equilateral 
type of non-Gaussianity, for which the non-Gaussian halo bias has only a very weak scale 
dependence. Note that if we find that the analytic modelling of the non-Gaussian bias works 
for relatively large non-Gaussianity, it will also be valid for small values of /nl- 



6 Results 

In this section we present the outcome of our N-body simulations of different types on non- 
Gaussianity. Here, we used the standard templates (namely the local, equilateral and or- 
thogonal templates) to create the initial conditions for the suite of N-body simulations. The 
templates are useful as toy models to test and calibrate the analytical predictions. If the 
analytical expressions can correctly reproduce the bias behaviour for shapes that have such 
different fe-dependence in the squeezed limit, this lends support to their validity. 

Although the simulations were targeted to study the non-Gaussian halo bias, we first 
consider the effect of non-Gaussianity on the halo mass function. Afterwards, we analyse in 
detail the scale-dependent halo bias induced by the different types of non-Gaussianity. 

6.1 Halo mass function 

Primordial non-Gaussianity affects significantly the number density of objects which corre- 
spond to the tails of the probability distribution function (PDF) of the initial density field 
[28] (for recent reviews, see [22, 23]). Galaxy cluster-sized halos originate from rare high 
peaks in the density field smoothed on the scale of the cluster mass. If the PDF has a posi- 
tive skewness on this mass scale, more of these high peaks are initially present and give rise 
to more galaxy clusters. 

The ratio of the non-Gaussian to the Gaussian halo mass function, i?^*^(M, z), can 
be modelled analytically in the framework of the Press-Schechter formalism [28, 29] and 
the excursion set approach [84, 85]. For the comparison with our N-body results, we use 
the formulas of [28] (MVJ) and [29] (LV), who used the saddle-point technique and the 
Edgeworth expansion, respectively, to derive the following expressions 

RmvAM.z) = exp . , ,o /o A^ + V ^ ~ Ac(z)53,m/3 

V 6(7^./ / 6^1 - Ac(z)53,M/3dlncJM V 

(6.1) 
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Figure 5. Fractional difference in the mass function derived from non-Gaussian simulations of the 
local type and from Gaussian simulations. Different realizations are depicted with symbols of different 
shapes. For clarity, the data points of different realizations are shifted slightly along the x-axis. The 
solid and dashed lines show the predictions of [28] (MVJ) and [29] (LV) using the best-fit fudge factor 
of g = 0.75 and q — 0.9, respectively. 



and 

6Ac{z) [ ' V 4f Ki J din aM \ aj^ J\ 

where aj^j is the linear mass variance on a mass scale M at z = 0. The redshift dependence 
is modelled by the collapse threshold, Ac{z) = ^6cD{z = 0)/D{z), which includes a fudge 
factor, y^^^. The non-Gaussianity is described by the normalized skewness of the linear den- 
sity field, 6, on a mass scale M at z = 0, S^^m = {^'Ii)/^m- The normalized skewness scales 
linearly in /nl- The magnitude and sign of S^^m depend on the shape of non-Gaussianity, 
while the mass dependence is only weakly dependent on the type of non-Gaussianity. 

In Fig. 5, we show the fractional difference in the non-Gaussian and Gaussian mass 
functions derived from the simulations of the local type of non-Gaussianity. The four panels 
correspond to different redshifts, z = 1.5 (top left), z = 1 (top right), z = 0.67 (bottom left), 

^^This fudge factor was introduced in [46] and interpreted as a modification to the collapse threshold which 
also applies to the Gaussian mass function. Here, however, we do not adopt this interpretation and consider 
q only as a fudge factor in the ratio of the non-Gaussian to the Gaussian mass function, which we calibrate 
using N-body simulations. 
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Figure 6. Effect on tfie halo mass function induced by the orthogonal shape of non-Gaussianity. 
Symbols and lines have the same meaning as in Fig. 5. Due to small signal to noise, the best-fit value 
for q at z = is only determined to about 20%. 



and z = (bottom right). The filled and open symbols represent the results for /nl = 250 
and /nl = 60, respectively. The corresponding analytic predictions of MVJ and LV (with 
the best-fit value for q) are depicted by the solid and dashed lines, respectively. Overall, 
the agreement between simulations and predictions is very good. Only at the low-mass 
end at z = 0, the data points do not follow the theory lines very closely. This is not of 
practical importance, as the difference is much smaller than observationally uncertainties 
in halo mass measurements. The applied fudge factor of q = 0.75 for MVJ and q = 0.9 
for LV is in agreement with our findings in [48]. Other groups [45, 46], who worked with 
friends-of- friends (FOF) halos, used q ~ 0.75 for both predictions, while [44], who applied 
a spherical-overdensity (SO) halo finder as we do, found g ~ 1 for LV and q « 0.75 for 
MVJ (V. Desjacques, private communication). For the standard definitions of FOF and 
SO halos, i.e. using a linking length of 0.2 times the mean particle distance (FOF) and 
the redshift-dependend virial overdensity (SO), [23] found that for a given halo mass the 
effect of non-Gaussianity on the number density of FOF halos is smaller than on the number 
density of SO halos. Note however that it remains to be seen whether FOF or SO halos most 
closely match observationally-selected halos (e.g., SZ-selected or X-ray selected clusters). 
Nevertheless, the fact that q is quite close to unity in both cases, suggests that this will not 
introduce a dominant systematic effect. 

The mass function results from our non-Gaussian N-body simulations of the orthogonal 
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Figure 7. Fractional difference in the mass function caused by non-Gaussianity of the equilateral 
type. The triangles correspond to two simulations with the same realization of the initial Gaussian 
field but differing in the method how the initial non-Gaussian contribution was computed (Eq. 4.5 
(black) and Eq. 4.9 (red)). The green squares depict the results of a second simulation with a different 
realization of the initial Gaussian field. Dashed and solid lines show the analytic predictions of [29] 
and [28]. 

type are presented in Fig. 6. Here, the /nl values are —1000 (filled symbols) and —250 
(open symbols). Note that these simulations used a rather small grid size {Ng = 400) for 
the non-Gaussian initial density field. Hence, halos with masses below 5 X IO^^Mq/Zi are 
expected to show a smaller non-Gaussian effect than predicted by the theory. This is indeed 
the case, as can be seen in the figure. More interestingly, even at the high-mass end the 
non-Gaussian effect is smaller than one could expect. To bring the theory in agreement with 
the N-body data, the fudge factor has to be quite small, decreasing from q 0.75 to q ^ 0.5 
with decreasing redshift. 

Lastly, we consider the increase in the halo mass function due to non-Gaussianity of the 
equilateral type. In Fig. 7, the results of three simulations with /nl = 1000 are shown. Two of 
them, "eql. a)", started from initial conditions computed with our modified ansatz, Eq. (4.9), 
using a non-Gaussian grid size of 400. The third simulation, "eql. b)" (black triangles), used 
initial conditions generated with our original method, Eq. (4.5), with Ng = 1024 (see also 
App. A). Remember that in the case of the equilateral type, artificial contributions to the 
power spectrum are kept under control also when our original ansatz is used. As expected, 
the N-body results agree with each other for halo masses above 5 x 10^^ Mq/H. The fall-off 
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of "eql. a)" at the low-mass end is caused by the smah non-Gaussian grid. Interestingly, we 
find that, in the case of the equilateral shape of non-Gaussianity, q is very close to one and 
does not change significantly with redshift. 

In conclusion, on the mass scales probed, the non-Gaussian effect on the halo mass 
function can reasonably well modelled by the analytic predictions, if a shape-dependent 
fudge factor is taken into account. In the case of the orthogonal type, there is a hint that q 
is in addition redshift dependent. 

6.2 Halo bias 

In order to explicitly explore the dependence on the halo mass, we divide the halos into five 
mass bins: 



3 X lO^-'^ Mrs/h < Mhaio < 6 X 10^^ Mq/H 

6 X lO^^Mo/Zi < Mhaio < 1-2 X IO^^'Mq/Zi 

1.2 X W^^Me/h < Mhaio < 2.4 x IO^^Mq/Zi 

2.4 X IO^^Mq/Zi < Mhaio < 4.8 x IO^^'Mq/Zi 

4.8 X lO^^Mo/Zi < Mhaio < 9.6 x IO^^Mq/Zi. 



The number of halos per mass bin decreases with mass. In particular, at high redshift and 
for large masses the exponential fall-off in the mass function reduces the number of halos per 
bin quickly. In the following, we only consider mass bins with more than 250 halos. 

In order to keep the statistical error as small as possible in our analysis, we compute 
and fit the effect of the non-Gaussian halo bias in the following way. First, we compute 
the difference between the non-Gaussian and Gaussian bias from the corresponding N-body 
simulations: 

^b{k) - p^^^^^ p^^^^ . (6.3) 

As the non-Gaussian and the Gaussian simulation share the same realization of the 
initial Gaussian field, A6(Zc) is almost free of sample variance and, in addition, has smaller 
shot noise than h^'^{k) and lfi{k) individually. We estimate the error on A6(A;) directly from 
the distribution of A6(k) in each k bin (for details, see App. B). 

Following Eq. 2.1, we model A6(Zc) by 

where is the linear halo bias obtained from the Gaussian simulation on large scales (see 
App. B). The scale- independent shift, A6/, can be predicted by the difference in the linear 
bias derived from the non-Gaussian and Gaussian mass functions using the peak-background 
split approach [44] 

Ode 

with R^^{M) being the ratio of the non-Gaussian to the Gaussian mass function. Here, 
however, we treat A6/ as a free parameter and compare it later on with the prediction 
derived from the mass functions. We choose q as the second free parameter. All other 
quantities in Eq. (6.4) are computed from the theory and are kept fixed. 
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Figure 8. Local non-Gaussian bias of lialos with mass 1.2 x 10^^Mpc//i < Af < 2.4 x 10^**Mpc//i 
at z 0. The difference of the bias measured from the non-Gaussian simulations and Gaussian 
simulations is depicted by the data points (for computation of the error bars, see App. B). The 
solid lines show the best fit using the model given in Eq. 6.4. The dot-dashed lines show the model 
predictions when the scale-independent bias shift, A6/, is neglected. The short-dashed lines neglect 
the term which is non-linear in /nl (see text for details). Thick lines and red symbols correspond to 
/nl = 250, while thin lines and blue symbols show the results for /nl = 60. Note that we actually 
show Ab + 0.1 to allow for a logarithmic scale. 

In Fig. 8, we show as an example the effect of local non-Gaussianity on the halo bias 
for halos of mass 1.2 - 2.4 x IO^^Mq/Zi at z = 0. Note that we plot Ab{k) + 0.1. As Abi is 
negative, this addition of 0.1 is needed to still make use of the logarithmic scale. 

The different line types visualize the effect of the different terms in Eq. (6.4). The solid 
lines show the best fit to the data (using all modes up to /cmax = 0.1 Mpc//i) and include all 
terms given above. The short-dashed lines neglect Abj appearing inside the square brackets 
in Eq. (6.4). The inclusion of this term makes the non-Gaussian bias non-linear in /nl 
[59, 61], since Ab[ depends on /nl- The dot-dashed hues neglect Abj completely. This 
scale-independent bias shift becomes important on smaller scales (k > 0.02), for which the 
scale-dependent part becomes small. 

An example of the measured non- Gaussian bias from the simulations of the orthogonal 
type is given in Fig. 9. Here, the halos have again a mass of 1.2 — 2.4 x 10^^ Mq//i, but were 
found in simulation outputs at z = 1. Consequently, the number of halos is smaller than 
in the previous figure and the residual shot noise is larger. The line types have the same 
meaning as before. On large scales, the halo bias scales as ~ as predicted by the theory. 
Hence, with increasing wavenumber, the effect does not drop as rapidly as in the local case 
and extends to smaller scales. 

Next, we present an example of the non-Gaussian halo bias induced by the equilateral 
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10 




Figure 9. Same as in Fig. 8, but for the orthogonal shape of non-Gaussianity. Note, however, that 
here the redshift of the halos with mass 1.2 x lO^^Mpc/h < M < 2.4 x 10"Mpc//i is z = 1. 
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Figure 10. Same as in Fig. 8, but for the equilateral shape of non-Gaussianity. Note the linear scale 
of the y-axis. 



shape. In Fig. 10, the halo bias corresponding to two different mass bins at z = is shown. 
As expected, the scale dependence is very weak and in agreement with the theoretical predic- 
tions. In particular, the observed mass dependence of the effect is consistent with the model 
predictions. 





Figure 11. Best-fit values (fitting all modes with k < 0.1Mpc//i) of the fudge factor q (left) and 
the scale-independent shift, A6/, normalized by /nl (right). Different colours correspond to different 
redshifts (see key in right panels) . The symbol coding is the same we used in the mass function plots 
(Fig. 5, 6, and 7), see text for details. On the left-hand side, all points in between two vertical dashed 
lines correspond to the same halo mass bin, but show different redshifts (coulor) , realizations (symbol 
shape), and /nl values (open/filled symbol). The horizontal dashed lines show the q values estimated 
from the LV non-Gaussian mass function. The coloured solid lines on the right-hand side show the 
predictions for the scale-independent bias shift using Eq. 6.5. 
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After having discussed for each type of non-Gaussianity typical examples, we show the 
complete set of best-fit values of the fitting parameters in Fig. 11. On the left-hand side, the 
best-fit values of the fudge factor q are presented. Different colours correspond to different 
redshifts: z = 1.5 (red), z = 1 (green), z = 0.67 (blue), and z = (magenta). Triangles, 
boxes, and circles depict the three different realizations of the initial Gaussian random field 
used for the generation of the initial conditions. In the case of the local and orthogonal type, 
the open symbols correspond to /nl = 60 and /nl = —250, respectively. Filled symbols show 
the results for /nl = 250 (local), /nl = —1000 (orth.), and /nl = 1000 (eql.). For clarity, 
the points of each mass bin are spread over the range of each mass bin. 

For the local type of non-Gaussianity, we recover within the error bars the same q value, 
which was needed to bring the mass functions in agreement with analytic predictions. For 
the orthogonal shape, the effect is suppressed by a factor which is similar to the one we 
found for the non-Gaussian mass function, q ~ 0.6. This suppression of the non-Gaussian 
halo bias effect is clearly redshift and halo mass dependent. In the case of the equilateral 
type of non-Gaussianity, the error bars are as expected very large. Nevertheless, also in this 
case there is a hint of a redshift and halo mass dependence of the effect. Note, however, that 
this suppression visible for the non-Gaussian halo bias effect for the equilateral shape was 
not found for the effect on the mass function, for which the measured fudge factor was g ~ 1. 
These findings are very interesting and — if solidified by larger simulations — may help to 
lead to a better theoretical understanding of the halo biasing (see discussion in [72]). 

On the right-hand side of Fig. 11, the best-fit values of Abj normalized by /nl are 
shown. The colour and symbol coding is the same as before. As open and filled symbols 
(corresponding to different /nl values) are consistent with each other, we can infer that the 
scale-independent shift is linear in /nl for the /nl values probed. The solid lines represent 
the predictions from Eq. (6.5) using the LV mass function ratio, R^{M)^ and taking the 
measured fudge factor q from the mass function into account. Keeping in mind that the LV 
mass function is not in all cases a good fit to the mass functions derived from our simulations 
(see for example bottom right panel of Fig. 6) , the agreement between the best-fit values and 
the predicted bias shift is very reasonable. 

7 Conclusions 

N-body simulations are an indispensable tool to test, develop and calibrate any statistical 
analysis of large-scale structure. In this paper, we have further explored the issue of setting 
up generic non-Gaussian initial conditions for N-body simulations. In a previous paper [48] 
we began addressing this issue, focussing — as we do here — on non-Gaussianities specified by 
a non-zero bispectrum. In [48] we focussed on the non-Gaussian corrections to the halo mass 
function and to the non-linear matter power spectrum, which, for the scales and redshfits of 
interest, could be reliably computed from simulation boxes of 600Mpc//i on a side. Here we 
concentrate on the large-scale non-Gaussian halo bias. The scale-dependent bias is evident 
on very large scales thus needs multiple simulations of much larger boxes to be measured 
reliably for /nl values not too large as to push beyond the regime of validity of the analytical 
description of the effect. 

The scale-dependent non-Gaussian halo bias has been recognized to be a very com- 
petitive approach to constrain primordial non-Gaussianity (e.g., [37, 40, 41, 59]) yielding 
forecasted error bars on the local non-Gaussian parameter of A/^^ ~ 1 which makes this 
approach formally better than an ideal CMB experiment [86]. So far this effect has been 
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explored with N-body simulations only for the so-called local type of non-Gaussianity. How- 
ever, the shape of non-Gaussianity is a window into the physical mechanism of inflation and 
the generation of primordial perturbations, and it is of value to be able to simulate initial 
conditions for general bispectrum shapes. 

The technique presented in [48] was not suitable to study the halo bias because, in 
general, higher-order contributions would leave an artificial signature on the large-scale power 
spectrum (which are the scales where the scale dependence of the halo bias is evident). These 
components could be kept under control only in specific cases. Here the issue is solved 
and these components are always under control. This however comes at the expense of 
computational cost: the new expressions are not separable (even if the bispectrum itself is) 
meaning that computational short-cuts that could be implemented for separable bispectra 
before, now cannot be applied in the same way. 

Nevertheless we find a speed-up solution which involves sampling only the Gaussian 
contribution to the potential with a full resolution grid and using a smaller grid for the 
non-Gaussian part. We assess the performance of this approach for the local case, where 
the computation can be performed in real space and is not computationally intensive. By 
using a smaller grid size for the non-Gaussian part of the initial conditions, we can keep the 
computational time taken by the calculation of the initial condition comparable to the time 
needed to run the simulation. 

The form of the bispectrum predicted by inflationary models can be complicated and 
in general is not-factorizable. As factorization is very important for an efficient analysis of 
CMB data, physical bispectrum shapes have been approximated by factorizable templates. 
We highlight here that, since the halo bias is dominated by squeezed configurations, what is 
a good template for CMB analysis may not be correct for the halo bias. We show that the 
equilateral template results in roughly the correct scale dependence of the halo bias for the 
corresponding physical models, but with a shift in amplitude; if unaccounted for, it could 
introduce a bias on the estimated /nl of up to 50%. Fortunately, this shift can be quantified 
well. Enfolded and orthogonal templates on the other hand do not work as well. They both 
lead to a non-Guassian bias that scales as on large scales. This is different from the halo 
bias induced by the corresponding physical models: The scale dependence of the halo bias 
caused by modified-initial-state type of non-Gaussianity is almost degenerate with the one 
caused by the local type of non-Gaussianity. The halo bias effect of the orthogonal physical 
shape is very similar to the effect of the equilateral shape. 

We performed several simulations with non-Gaussian initial conditions of different shapes 
to further test and calibrate the analytic predictions. For this purpose we used the standard 
templates: these are useful toy models as they span a range of scalings in the squeezed regime. 
If the analytical expressions can correctly reproduce the bias behaviour for shapes that have 
such different /c-dependence in the squeezed limit, this lends support to their validity. 

Revisiting the mass function predictions, we find that the non-Gaussian effect on the 
halo mass function can be well modelled by the analytic predictions, if a shape-dependent 
fudge factor, q, is taken into account. In the case of the orthogonal type, there is a hint that 
this factor depends on redshift and halo mass. 

For the halo bias we also find that the analytic predictions work well once a a shape- 
dependent normalization factor q is included: q effectively re-scales /nl- Interestingly, for 
the non-local types of non-Gaussianity, this fudge factor varies with redshift and halo mass. 

Throughout we pay particular attention to the operational definition of non-Gaussian 
bias, taking into account that a non-zero /nl modifies the mass function and thus modifies 
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also the linear bias. However, a Gaussian distribution with the same mass function would 
have a linear halo bias modified in the same way. 

In most inflationary models considered so far, the bispectrum scales either as the equi- 
lateral shape (~ or as the local shape (~ k~^) in the squeezed limit (see, however, the 
quasi-single field models of [87, 88], which have intermediate scalings). The equilateral shape 
leads to a small and almost scale-independent halo bias and is thus not easily accessible 
with this technique. Other shapes that, on large scales, have the same scale dependence as 
the local one, may in principle be distinguished through the halo bias dependence on halo 
mass on intermediate scales. On large scales, however, the limits for the local /nl can be re- 
interpreted as limits on e.g., modified-initial-state type of non-Gaussianity, by an appropriate 
rescahng of /nl: f^£^'^- ~ Sf^^. 
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A Simulations of the equilateral type using two different approaches 

In this section we compare the results of two N-body simulations with /nl = 1000 non- 
Gaussian initial conditions of the equilateral type. Both simulations have a box size of 
1875Mpc//i and a particle load of 1024^. The simulations were stB-rtcd. dit ^initial 

= 49 and 

carried out with Gadget-2 using the same numerical settings. In both cases, halos were iden- 
tified with AHF [82] . The only difference between the simulations is the different prescription 
and implementation how the non-Gaussian part of the initial gravitational potential was gen- 
erated. In the first case, "equilateral a)", we used our modified ansatz for Eq. (4.9), 
and a grid size of 400 and 1024 for the non-Gaussian and Gaussian part, respectively. The 
initial conditions of the second simulation, "equilateral b)", were generated using our original 
ansatz, Eq. (4.5). The integral appearing in this ansatz is a sum of convolutions and can 
be computed efficiently with the help of FFTs. Hence, in this case, we were able to use an 
equal large grid size of 1024 for the non-Gaussian part, too. In both cases, we use the same 
realization of the Gaussian random field. Note although the two approaches are different at 
the level of they both produce a non-Gaussian field, which has, at the leading order, the 
same power spectrum and bispectrum. However, high-order spectra, e.g. the trispectrum, 
are different. 

In the remainder of this section, we consider observables which are affected by primor- 
dial non-Gaussianity, and analyse if the different ways of setting up the initial conditions lead 
to differences in these quantities. We start with the non-linear matter power spectrum. In 
Fig. 12, we show the ratio of the power spectrum measured from the two non-Gaussian simu- 
lations to the power spectrum obtained from a Gaussian simulation with the same realization. 
The cyan/blue and magenta/red lines correspond to redshift z = and z = 1, respectively. 
The vertical black line shows the Nyquist frequency of the 400^ grid used for the non-Gaussian 
part in "equilateral a)" , which corresponds to the smallest scale still resolved by the grid. For 
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Figure 12. The ratio of tlie non-linear power spectrum obtained from the two non-Gaussian simu- 
lations to the Gaussian power spectrum at redshift z = 1 and z = 0. The black vertical line shows 
the Nyquist frequency, up to which the primordial non-Gaussian field is sampled in "equilateral 
a)". 

both redshifts, significant differences between the non-Gaussian simulations only appear on 
scales smaller than the smallest scale sampled by the non-Gaussian grid used in "equilateral 
a)", i.e. on the right hand side of the vertical line. On the largest scales, the small number 
of modes leads to differences in the power spectra, since the term ($'^$^'^) oc ($'^$'^<I)'^) 
does not completely vanish. The very good agreement on the intermediate scales confirms 
the theoretical expectation, that the non-Gaussian amplification of the power spectrum is at 
leading order caused by the primordial bispectrum. 

Next, we show in Fig. 13 the halo bias for different halo masses at redshift z = 0. The 
halo bias is computed from the ratio of the halo-matter cross power spectrum to the matter 
power spectrum, see Eq. (5.2). The dashed and solid lines correspond to the simulations 
"equilateral a)" and "equilateral b)", respectively. As the effect of non-Gaussianity of the 
equilateral type on the halo bias is very weak (see Fig. 2), no effect of spatially unresolved 
non-Gaussianity, like we observed in the local case (see Fig. 3), is visible here. The small 
differences between the two simulations can be ascribed to the differences in the shot noise. 

Lastly, we consider the halo mass function derived from the two non-Gaussian simu- 
lations. The ratio of the halo number densities measured from the simulations "equilateral 
a)" and "equilateral b)" as a function of halo mass is shown in Fig. 14. Similar to our test 
with the local type of non-Gaussianity (see Fig. 4), the number of halos with mass below 
3 X 10^^ Mq/Zi is reduced by a few percent in the simulation for which a smaller non-Gaussian 
grid was used to set up the initial conditions. Halos with masses larger than ~ 5 x 10^^ Mq/H 
are not affected by the smaller grid size. Both, the magnitude of the effect and the affected 
mass range are a bit larger than in the local case, which is probably due to the larger effect 
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Figure 13. Halo bias for different halo masses in units of Mq//i at redshift z — 0. The dashed 
and solid lines show the results obtained from the non-Gaussian simulations "equilateral a)" and 
"equilateral b)" , respectively. 



on the halo mass function caused by the primordial non-Gaussianity with = 1000 instead 
of f^l = 250 (see Fig. 5 and Fig. 7). 

As in the case of the non-linear matter power spectrum, no additional effects due to 
the differences in the initial non- Gaussian fields of "equilateral a)" and "equilateral b)" are 
noticeable. This strengthens the analytic modelling that the effects on the mass function due 
to non-Gaussianity are primarily caused by the primordial bispectrum. 



B Details on the fitting of the non-Gaussian bias 

We assume that the modes of the halo density field are related to the matter density by 
(5h = 6(A;)(5m -|- n{k), where n{k) is a noise variable which models the stochasticity in the 
halo formation process and the discrete nature of halos. The quantity which we want to 
measure from the simulations is the bias b{k). Assuming that the noise is uncorrelated with 
the density the bias is given by 

In practice, the average is taken over the finite number of modes per k bin. The term 
(Re(n(5^)) is therefore not exactly zero. Assuming Re{n6^) has a Gaussian distribution with 
variance cr^^, the error of b{k) is then given by / V -^modcs / (I (^m P ) , where A^modes is the 
number of independent modes. 

In Fig. 15, we show the square root of the variance of the noise distribution for halos at 
z = 1 with masses of 1.2x 10^^ — 2.4x IO^^Mq/Zi, measured from different simulations. Filled 
triangles and open squares depict the estimates derived from Re(n5j^)/|(5m| and lm{n5^) /\Sra\, 
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Figure 14. The ratio of mass function of halos found in the simulation "equilateral a)" and in the 
simulation "euilateral b)" at redshift z = 1 and z = 0. The error bars show the Poisson error. For 
clarity, the green circles are slightly shifted along the x-axis. 

respectively. The dashed lines show the Poisson noise prediction cipoisson = -\/l/(2nh), where 
rih is the halo number density. For the halos shown, the measured noise is in good agreement 
with the Poisson noise prediction. However, in general, depending on the halo mass and 
redshift, the actual noise can be smaller or larger than the Poisson noise and can also be 
scale-dependent (e.g., [89]). 

As the Gaussian and non-Gaussian simulations share the same realization of the Gaus- 
sian field, we expect that in the quantity 

, _ Re«^^^^*) Re{6^6g*) 

the corresponding noise terms Re{n^^5^'^*)/\6^'^\ and Re{n'^6^*)/\6^\ cancel each other 
partly. This is indeed what we find in the simulation data. The blue symbols in Fig. 15 
represent cr„ derived from the distribution of A in each fc-bin. The solid line is a linear fit 
to the data. When we fit the non-Gaussian halo bias models to the simulation data, we use 
this fit as our noise estimate. The resulting x^/d.o.f. values are close to one. 

Another ingredient in our fitting procedure of the non-Gaussian halo bias is the Gaussian 
linear bias, bf. We estimate the linear bias by fitting the bias derived from the Gaussian 
simulations on large scales. In Fig. 16, the bias measurements for halos of different masses 
at redshift z = 1 are shown. The dashed lines depict the corresponding best fits of the linear 
bias. Only data points on large scales, k < 0.07, were included in the fitting. 

The remaining quantities needed for the fitting of the non-Gaussian bias (see Eq. (6.4)), 
like the growth function, the form factor, etc., are computed numerically for the given redshift. 
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Figure 15. Variance of tlie distribution of Re{n6^)/\S^\ and lm{nS^)/\Sm\ derived from N-body 
simulations for halos with mass of 1.2 x 10^^ — 2.4 x IO^^Mq/H at 2 = 1. Red and green symbols 
correspond to a Gaussian and non-Gaussian simulation of the local type with /nl = 250, respectively. 
The blue symbols show the noise derived from the distribution of A (see. Eq. B.2). The dashed 
horizontal lines depict the Poisson noise prediction. The solid line shows the linear fit to the data. 

type of non-Gaussianity, and halo mass^'^. 

Having all quantities at hand, we fit the non-Gaussian bias derived from the simulation 
data with the model given in Eq. (6.4) out to wavenumbers k < k^^x- We use two different 
fitting ranges, a conservative /cmax = 0.03Mpc//i and 0.1Mpc//i. 
Both fitting ranges yield results consistent with each other. However, in the case of the 
equilateral and orthogonal shape, the errors on the fitting parameters are significantly smaller 
(approximately by a factor of ten and two, respectively) when the larger fitting range is 
used. For the local type of non-Gaussianity, the larger fitting range reduces the error by 
approximately 30%. 
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